Hydrologic Design 4501
Chapter 7: Unit Hydrograpgh and Routing
Ardeshir Ebtehaj
University of Minnesota
Table of Contents
Chapter 7: Unit Hydrograpgh and Routing
1- Linear Reservoir and Convolution
1-1 Concept of Convolution
1-2 Instantaneous Unit Hydrograph (IUH)
1-3 Unit Hydrograph
2- Computation of the Unit Hydrograph
2-1 Forward substitution and least-squares
2-2 S-Hydrograph
2-3 Synthetic Unit Hydrographs
2-3-1 Concentration time
2-3-2 Time to Peak
2-3-3 Peak Discharge
2-3-4 Construction of SCS synthetic UH
4- Flow Routing
4-1 Level Pool Routing
4-2 Hydrologic River Routing
1- Linear Reservoir and Convolution
A linear reservoir is a simple and meaningful representation of the watershed response to precipitation events. Continuity equation for a linear reservoir is written as follows:
Figure 1: Schematic of a linear reservoir. A reservoir is linear if we assumed that the storage function 
is linearly related to the outflow through the constant K.
is linearly related to the outflow through the constant K.where 
is the outlfow, 
is the inflow 
(rainfall excess), 
is change in the basin storage over time, and K represents a time scale 
that captures the residence time of the reservoir. Here, the storage is the sum the rainfall excess minus the outflow up to time t.
is the outlfow, is the inflow (rainfall excess), is change in the basin storage over time, and K represents a time scale that captures the residence time of the reservoir. Here, the storage is the sum the rainfall excess minus the outflow up to time t. Now upon substitution, we obtain and have to solve a first-order ordinary differential equation (ODE):
Multiplying both sides by 
, one can get:
, one can get:
Now because of the product rule, we can rewrite the above equation as follows:
,and thus, we have,
Therefore, we can obtain the outflow of a linear reservior as follows:
1-1 Concept of Convolution
The above solution states that the outflow at any time is a linear sum of the contributions by an initial flow (first term) that goes to zero over time and the integral term (second term) that represents the effects of the reservoir resident time Kon the outflow. The integral term is a convolution integral. Learning the concept of convolution is crucial for understanding the concept of the unit hydrograph.
Mathematical convolution for two generic functions 
and 
is defined as follows:
and is defined as follows:
Figure 2: Visualization of two examples of the convolution integral. We can see that the convolution of a function with a kernel 
may result in translation (shift) and dilution of . We will show later on the the rainfall excess can be represented by and the kernal function or the unit hydrograph is a characteristic of the watershed. The result of the convolution of excess rainfall with the unit hydrograph will the be watershed outflow.
with a kernel may result in translation (shift) and dilution of . We will show later on the the rainfall excess can be represented by and the kernal function or the unit hydrograph is a characteristic of the watershed. The result of the convolution of excess rainfall with the unit hydrograph will the be watershed outflow. 1-2 Instantaneous Unit Hydrograph (IUH)
Now looking back on the integral in the linear reservoir outflow equation, we can see that the inflow 
resembles 
and the convolution kernal 
resembles 
resembles and the convolution kernal resembles 
and 
. If we define the input rainfall excess as a unit impulse function as follows:
Then, we can define the impulse response function of a linear reservoir as follows:
The function is often called the instantaneous unit hydrograph (IUH), which is the response of a linear reservoir to a unit excess rainfall with an infinitesimally small duration.
is often called the instantaneous unit hydrograph (IUH), which is the response of a linear reservoir to a unit excess rainfall with an infinitesimally small duration.1-3 Unit Hydrograph
Similar to the above formalism for a linear reservoir, a watershed has a response function to a unit pulse of excess precipitation --- albeit a more complex one. This response function is known as the unit hydrograph, which is a time-invariant characteristic of a watershed. More specifically, the unit hydrograph 
of a watershed is defined as tthe direct runoff hydrograph (DRH) resulting from a unit pulse of excess rainfall, which is defined as a unit depth of excess rainfall (i.e., 1 
, 1 
, etc) distributed uniformly over the entire watershed for an effective duration of 
.
of a watershed is defined as tthe direct runoff hydrograph (DRH) resulting from a unit pulse of excess rainfall, which is defined as a unit depth of excess rainfall (i.e., 1 , 1 , etc) distributed uniformly over the entire watershed for an effective duration of .
Figure 3: The response of a linear reservoir to a unit impulse (left) resembles the response of a watershed, called unit hydrograph, to a unit pulse of rainfall excess (right), which is a unit amount of rainfall, uniformly distributed over the watershed for duration of .
.Note that convolution is a linear operator
which means that we can superimpose the impacts of two pulses of excess rainfall over a watershed through linear combination of their individual outflow responses.
Figure 4: Superimposition of the response functions (i.e., re-scaled and shifted unit hydrographs) to excess rainfall hyetographs. For example, the above schematic shows the response to two pulses of excess rainfall, where 
.
. Therefore, having the unit hydrograph, one can obtain the outflow of a watershed as convolution of excess rainfall 
with the unit hydrograph as follows:
with the unit hydrograph as follows:
Below is an example of hydrograph superposition used to obtain the outflow of a watershed in response to a 24-hr storm -- when a 6-hourly unit hydrograph is available.
Figure 5: Example of 6-hr Unit Hydrograph convolution (Credit: The COMET Program).
2- Computation of the Unit Hydrograph
Figure 6: Convolution of excess rainfall hyetograph (ERH) at 
intervals with the unit hydrograph of a basin, sampled at 
points. The convolution results in a Direct Runoff hydrograph (DRH) at 
points, where 
.
intervals with the unit hydrograph of a basin, sampled at points. The convolution results in a Direct Runoff hydrograph (DRH) at points, where .Since, hydrologic data is not continuous, we must use a discretized version of the convolution equation as follows:
For the example shown in Figure 6, we have
where 
, 
. Therefore, the above convolution can be expressed as a linear system of equation as follows:
, . Therefore, the above convolution can be expressed as a linear system of equation as follows:
Through convolution, one can obtain the above set of equations through flipping the rainfall excess hyetograph and slide it over the coordinates of the unit hydrograph and then compute the sum over multiplications.
2-1 Forward substitution and least-squares
We can solve the above set of equations through forward substitution. In other words, one can obtain 
and then use the second equation to obtain 
and so forth.
and then use the second equation to obtain and so forth. One can arrange all data of the above linear system of equations in a matrix form as follows:
and solve it for 
values to obtain the unit hydrograph from the observed direct runoff hydrograph of a precipitation excess event.
values to obtain the unit hydrograph from the observed direct runoff hydrograph of a precipitation excess event. We can see that in the above linear system of equations we have more equations that unknown. Therefore, calculation of the unit hydrograph is often an overdetermined problem.
The above system of linear equation, which is resulted from a discrete convolution operation can be rewritten in matrix form as follows:
where 
is an 
column vector, 
is an 
matrix, and 
is an 
column vector.
is an column vector, is an matrix, and is an column vector. For the provided example, the matrix form is as follows:
The matrix form can be easily implemented into a software library for numerical linear algebra. As previously explained, the equation is an over-determined linear system, since we have more equations than unknowns. Two of the main methods used to solve this system of equations for are:
is an over-determined linear system, since we have more equations than unknowns. Two of the main methods used to solve this system of equations for are: 1. Forward/backward substitution: Write out the first or last L equations. Solve the first or last equations and sequentially substitute and solve the remaining equations. In this formalism only L equations are used, which result in multiple solutions to the problem.
2. Least squares: Using linear algebra, we can derive a formula to calculate the least squares estimate of the unit hydrograph ordinates as follows:
where 
is the matrix transposition operator, which flips the elements of a matrix over its diagonal elements. Note that by least squares, we mean that the solution has minimum error variance amoung all other possible solutions and the solution is unique.
is the matrix transposition operator, which flips the elements of a matrix over its diagonal elements. Note that by least squares, we mean that the solution has minimum error variance amoung all other possible solutions and the solution is unique. The above matrix formulation of the linear system can be solved efficiently using a technical programing language such as MATLAB. Note that, in the above formulation, we multiply both sides with 
to obtain a square matrix 
, which might be invertible. We need to note that is often a rectangular matrix with more rows than columns (
) and cannot be inverted.
to obtain a square matrix , which might be invertible. We need to note that is often a rectangular matrix with more rows than columns () and cannot be inverted. The matrix 
is called pseudo inverse of . In MATLAB, computation of the unit hydrograph can be simply implelmented using either of the following commands: u = P\q or u = pinv(P)*q.
is called pseudo inverse of . In MATLAB, computation of the unit hydrograph can be simply implelmented using either of the following commands: u = P\q or u = pinv(P)*q.--------------------------------------------------------------------------------------------
Example Problem 7.1: Find the half-hour unit hydrograph using the excess rainfall hyetograph and direct runoff hydrograoh shown in the following table.
Solution: Given 
, 
and since thus 
.
, and since thus .
% Least squares solution for Unit Hydrograph
clear
format short g
close all
t = 0:0.5:5;
% Direct runoff hydrograph (DRH)
q = [423 1923 5297 9131 10625 7834 3921 1846 1402 830 313]';
% Excess rainfall Hyetograph (ERH)
p = [1.06 1.93 1.81 0 0 0 0 0 0 0 0];
r = [1.06 0 0 0 0 0 0 0 0];
P = toeplitz(p,r);
% Unit Hydrograph
% u = inv(P'*P)*P'*q;
% or
u = pinv(P)*q;
%% Plots
yyaxis left
plot(t',q,'DisplayName','DRH','LineWidth',2)
hold on
plot(t(1:length(u)),u,'DisplayName','UH','LineWidth',2)
xlabel('\bf Time (hr)')
ylabel('\bf Discharge (cfs)')
yyaxis right
bar(t,p,'FaceAlpha',0.3,'DisplayName','Precipitation')
ylabel('\bf Precipitation excess(in)')
legend('show')
--------------------------------------------------------------------------------------------
Example Problem 7.2: Calculate the streamflow hydrograoh for a storm of 6 inch excess rainfall over three hours with hourly rate of p = [2,3,1] using the half hour unit hydrograph computed in the previous example and assume the base below is constant at 500 [cfs] throughout the flood. Check that the total depth of direct runoff is equal to the the total depth of excess precipitation for the watershed area of 7.03 miles square.
Solution:
%computation of the strearmflow hydrograph
clear
format short g
close all
dt = 0.5; % [hr]
t = 0:dt:5;
A = 7.03*5280^2; %[ft^2]
Q_base = 500; % [cfs]
u = [404, 1079,2343,2506,1460,453,381,274,173]; % [cfs/in]
p = [2,3,1]; % [in]
Q_direct = conv([0 0 u 0 0],p,'valid'); % [cfs]
Q = Q_direct+Q_base; % [cfs]
V_d = sum(Q_direct)*dt; % [cfs.hr]
V_d = V_d*3600/1; % [ft^3]
r_d = V_d/A*12; % [in]
Q_direct_1 = conv([0 0 u 0 0],[2 0 0],'valid');
Q_direct_2 = conv([0 0 u 0 0],[0 3 0],'valid');
Q_direct_3 = conv([0 0 u 0 0],[0 0 1],'valid');
yyaxis left
plot(t,repmat(Q_base,size(t)),'-ok','DisplayName','Base Flow','LineWidth',2)
hold on
plot(t,Q_base+Q_direct_1,'-or','DisplayName','Response to p = [2,0,0] [in]','LineWidth',2);
plot(t,Q_base+Q_direct_1+Q_direct_2,'-og','DisplayName','Response to p = [2,3,0] [in]','LineWidth',2);
plot(t,Q_base+Q_direct_1+Q_direct_2+Q_direct_3,'-ob','DisplayName','Response to p = [2,3,1] [in]','LineWidth',2);
xlabel('\bf Time (hr)')
ylabel('\bf Discharge (cfs)')
yyaxis right
bar(t,[p, zeros(1,8)],'FaceAlpha',0.3,'DisplayName','Precipitation')
ylabel('\bf Precipitation excess (in)')
legend('show')
2-2 S-Hydrograph
As we mentioned, the UH is a time invariant characteristic of a watershed in response to a unit pulse of excess rainfall with a specific duration 
. However, since the UH is obtained through a linear process, we can change a unit hydrograph from one duration to another. To that end, we need to derive the S-hydrograph based on the principle of superposition.
. However, since the UH is obtained through a linear process, we can change a unit hydrograph from one duration to another. To that end, we need to derive the S-hydrograph based on the principle of superposition. By definition, the S-hydrograph is the direct runoff hydrograph of a watershed resulting from a continuous excess rainfall at a constant rate of 1 [cm/hr] or [in/hr] for an infinite period of time.
In other words, the S-hydrograph 
is the response of the watershed to summation of an infinite number of unit pulse of excess rainfall with intensity 
and duration -- or a unit rainfall excess over an infinite period of time (Figure 7).
is the response of the watershed to summation of an infinite number of unit pulse of excess rainfall with intensity and duration -- or a unit rainfall excess over an infinite period of time (Figure 7). 
Figure 7: Conceptual representation of the S-hydrograph.
As a result, the unit hydrograph for 
duration, with intiger n values, can be obtained as follows:
duration, with intiger n values, can be obtained as follows:
Figure 8: Derivation of the unit hydrograph with duration 
, where 
, from an S-hydrograph that is generated from a unit hydrograph with duration .
, where , from an S-hydrograph that is generated from a unit hydrograph with duration .--------------------------------------------------------------------------------------------
Example Problem 7.3: Use 0.5-hr unit hydrograph in the following table to calculate 1.5-hr unit hydrograph.
Solution:
A shifted version of the hydrogrsph is shown between columns 2 and 3. Column 3 represents the sum of the unit hydrographs multiplied by 
[hr]. Column 4 shows the lagged unit hydrograph for 
[hr]
[hr]. Column 4 shows the lagged unit hydrograph for [hr]An excel file solving an example to derive a 1.5-h unit hydrograph from 0.5-hr unit hydrograph is posted in the course website. Column 5 computes the unit hydrograph for 
[hr].
[hr].% computing 1.5-hr hydrograph from 0.5-hr hydrograph using S-hydrograph
format short g
clear
dt = 0.5; % [hr]
dt_prime = 1.5; % [hr]
t=0.5:dt:6.0;
u_0_5 = [404 1079 2343 2506 1460 453 381 274 173 0 0 0];
M = zeros(length(t),length(t));
for i=1:length(t)
u_temp = u_0_5(1:end-i+1);
M(i:12,i) = u_temp(1:end)';
end
s_hydo = dt*sum(M,2);
s_hydro_lag = [0;0;0;s_hydo(1:end-3)];
u_1_5 = 1/dt_prime*(s_hydo-s_hydro_lag);
2-3 Synthetic Unit Hydrographs
Thus far, the methods we have described only apply for the basin with gage measurements of the streamflow, however, many basins are ungauged. Therefore, engineers have analyzed many unit hydrographs from gauged basins to obtain a universal functional representation for the unit hydrograph that can be generalized to ungauges basins. This representation of the unit hydrograph is often called Synthetic Unit Hydrograph (SUH) that can be applied to ungauged basins. These synthetic methods typically scale a general UH shape or calculate key ordinates based on easily determined characteristics of watershed and excess precipitation.
There are many of these methods, such as Snyder's or Clark's SUH; however, we will cover the very popular NRCS (SCS) dimensionless synthetic hydrograph method. The dimensionless SCS hydrograph and its triangular representation are shown in Figure 9 and tabulated in Table shown as Figure 11.
Figure 9: Dimensionless curvilinear and equivalent triangular SCS unit hydrographs (USDA, 1986)
Input parameters:
: duration of excess rainfall
: time lag between the centroid of rainfall excess hyetograph to the peak discharge
: time to concentration
: time to peak
: recession time
: base time for total duration of the scaled SUH
: peak flow of the scaled SUH.
2-3-1 Concentration time
As we previously discussed the time of concentration () for a watershed is the time for a raindrop to travel from the farthest point in the watershed to the outlet. The SCS recommends two methods for calculation of :
) for a watershed is the time for a raindrop to travel from the farthest point in the watershed to the outlet. The SCS recommends two methods for calculation of : - The lag method: defines the time lag as the time in hours from the center of the mass of the rainfall excess to the peak discharge, which can be computed as follows:
where 
is the average slope of the watershed in %, L denotes the hydraulic length of the main channel in feet, and 
[inches] represents the potential maximum retention that is obtained from the curve numbers (CN). The SCS method suggests the following relationship for the concentration time
is the average slope of the watershed in %, L denotes the hydraulic length of the main channel in feet, and [inches] represents the potential maximum retention that is obtained from the curve numbers (CN). The SCS method suggests the following relationship for the concentration time
- The upland or velocity method defines the concentration time as follows:
where 
and V are hydraulic length and velocity in [ft] and [ft/s], respectively. The velocity of the upland method can be obtained from Figure 8 based on the average slope of the watershed.
and V are hydraulic length and velocity in [ft] and [ft/s], respectively. The velocity of the upland method can be obtained from Figure 8 based on the average slope of the watershed. Clearly, if we discretize the the stream into k segments, the concentration time is the sum of travel time for different stream segments as we discussed before,
.
Figure 10: Velocities for upland method for estimation of (USDA, 1986).
(USDA, 1986).2-3-2 Time to Peak
Time to peak () is the time from the beginning of the rainfall to the time of the peak discharge, which can be expressed as a function of the time lag as follows:
) is the time from the beginning of the rainfall to the time of the peak discharge, which can be expressed as a function of the time lag as follows:
where is in 
, is the duration of the rainfall excess in , and is the lag time in .
is in , is the duration of the rainfall excess in , and is the lag time in . When is not known, the SCS method suggests
is not known, the SCS method suggests
and because 
, one can obtain the following relationship for approximating the time to peak:
, one can obtain the following relationship for approximating the time to peak:
2-3-3 Peak Discharge
From continuity equation, we know that the area under the unit hydrograph is equal to the volume of the excess unit rainfall - example 1 [in] all over the watershed. Here, the unit of in SCS method is the unit of discharge per depth of rainfall excess such as 
or 
.
in SCS method is the unit of discharge per depth of rainfall excess such as or .Therefore, base on conservation of mass of the depth of the rainfall and the triangular shape of the proposed SCS unit hydrograph, we have
which can be rearranged as
For the SI unit, in the above equation the area should be in 
and time is in second to obtain the discharge in . To convert the above formula to be used for watershed area in square kilometers and hourly time scale, we need to to multiply the above equation by 
.
and time is in second to obtain the discharge in . To convert the above formula to be used for watershed area in square kilometers and hourly time scale, we need to to multiply the above equation by . Moreover, the SCR method suggests 
and thus
and thus 
where 
. The 
value in English unit is 483.4, where A is in square miles.
. The value in English unit is 483.4, where A is in square miles. After obtaining the and , the curvilinear coordinates of the SUH can be obtained from the table in the next slide.
and , the curvilinear coordinates of the SUH can be obtained from the table in the next slide.
Figure 11: Ratios for dimensionless SCS synthetic unit hydrograph and mass curve.
2-3-4 Construction of SCS synthetic UH
- Compute the time of concentration () using the lag method or the upland-velocity method.
- Compute the time to peak
and then the peak discharge 
. - Compute the the base time of the SUH as
, which is 
for the triangular hydrograph and 
for the curvilinear one and the recession time 
. - Compute the SUH ordinates and plot them. For the the triangular SUH, the ordinates can be directly computed while for the curvilinear one the dimensionless ratios
shall be used.
Figure 12: The SCS unit hydrographs (USDA, 1986)
---------------------------------------------------------------
Example 7.2: Construct a 10-minute unit hydrograph for a basin with an area of 3.0 
and concentration time 
[hr].
and concentration time [hr].
[min] = 0.166 [hr] (duration of excess rainfall) 
[hr] (lag time)
[hr] (time to peak)
(peak flow discharge)
Figure 13: The calculated SCS SUH.
4- Flow Routing
Flow routing is the procedure of determining the downstream hydrograph from hydrographs at one or more points upstream.
As water enters a channel, either from hillslopes or upstream channels, the flow velocity shifts the location of the peak flow and surface roughness and storages along the water path attenuate the magnitude of the peak. The flow routing is referred to a class of methods that calculates the redistribution (widening) and translation of upstream hydrograph as a function watershed flow characteristics.
Figure 14: Left: Schematic of redistribution and translation (Chow et al. 1988). Middle: The storage accumulation and release during the flow routing. Right: Conceptualization of channels as storage reservoirs (COMET Program)
There are two main types of routing methods:
1. Lumped (Hydrologic) Routing: Flow is calculated at a particular point as a function of time alone.
2. Distributed (Hydraulic) Routing: Flow is calculated as a function of time and space.
For this course, we focus on the simpler lumped modeling.
Within lumped modeling, there are multiple methods depending on what system we are trying to model. We are going to focus on two common methods for flow routing in reservoirs (level-pool method) and stream channels (Muskingum method).
Derivation of routing methods begins with the continuity equation:
In routing,
is the known inflow hydrograph, while and are unknown. Just as in our discussion of linear reservoirs, we need a storage-discharge relationship, 
.
is the known inflow hydrograph, while and are unknown. Just as in our discussion of linear reservoirs, we need a storage-discharge relationship, . There are two types of storage-discharge relationships:
- (1) invariable reservoirs: The invariable are those in which there is a one-to-one relationship between the outflow and the reservoir storage. This scenario is valid in short and wide reserviors.
- (2) variable reservoirs: In variable reservoirs, the relationship between the outflow and storage capacity is not unique and follows a loop when the reservoir is filled or emptied. This nonlinear relationship mostly occurs in long and narrow reservoirs where the water velocity is significant and produces backwater effects in the reservoir.
Figure 15: Invariable reservoir (left) that is wide and deep compared to its length. The flow velocity is low and the water surface is almost horizontal. A variable reservoir (right) is a long, narrow and shallow reservoir with high flow velocity such that the water surface profile is markedly curved due to the backwater effect.
Figure 17: Photos of a reservoir, dam and a pond with hydraulic structures as flow controls.
4-1 Level Pool Routing
Level pool routing typically applies to man-made hydraulic structures like reservoirs, dams and stormwater ponds that are common in many watersheds. The key assumption is a horizontal water surface and negligible velocity in the reservoir. This allows an invariable relationship between storage (water surface height) and discharge, because any change in storage results in a uniform change in height of water surface.
Beginning with the continuity equation 
, we have 
, one can discretize it as follows assuming that the changes in inflow and outflow are linear over a small timestep :
, we have , one can discretize it as follows assuming that the changes in inflow and outflow are linear over a small timestep :
Figure 17: Discretization of the reservoir continuity equation over a small time interval.
The inflow 
and 
are known at all times while initial storage 
and outflow 
are specified as initial conditions, so we can re-arrange the equation into knowns and unknowns as follows:
and are known at all times while initial storage and outflow are specified as initial conditions, so we can re-arrange the equation into knowns and unknowns as follows:
To use this equation, we need the elevation-storage function 
and elevation-discharge 
relationships of the reservoir to create a storage-outflow function that allows obtaining the outflow at time step 
from the left-hand side of the above equation 
iteratively.
and elevation-discharge relationships of the reservoir to create a storage-outflow function that allows obtaining the outflow at time step from the left-hand side of the above equation iteratively. - The elevation-storage function is simply constructed from the geometry of the reservoir or pond.
- The elevation-discharge is slightly more difficult to obtain; however, the outflow of storage facilities is often controlled by hydraulic structures such as sluice gates, weirs, pipes, and orifices which typically have well-defined equations based on thr water head (Figure 17).
For example, for an uncontrolled ogee crest spillway, the outflow as a function of total water head on the crest is:
is the effective length of the crest
is the discharge coefficient 
is discharge in [cfs].
Figure 18: Discharge relationship for vertical-faced ogee crest (US bureau of reclamation 1987). The figure is for English system and 
is the water surface upstream from weir drawdown.
is the water surface upstream from weir drawdown.For a general reservoir, once the storage-outflow function is calculated, as schematically shown in Figure 19, one can calculate the RHS of Eq.(1) and compute the LHS.
Figure 19: Change of storage during a routing period and development of the storage-outflow function for level pool routing on the basis of storage-elevation-outflow curves (Chow et al. 1988).
---------------------------------------------------------------
Example 7.3: Inflow hydrograph is specified at 
min. Horizontal area of the reservoir is 
acre (43,560 
) and 
. Given the inflow hydrograph, compute the outflow hydrograph.
min. Horizontal area of the reservoir is acre (43,560 ) and . Given the inflow hydrograph, compute the outflow hydrograph.
Figure 20: The left four columns represent the problem inputs and construction of the storage-outflow curve 
. The last six columns represent the inflow and the routed outflow. An excel worksheet is posted in the class website.
. The last six columns represent the inflow and the routed outflow. An excel worksheet is posted in the class website.Solution:
- First we need to construct Q (y-axis, col 2) versus
(x-axis) (col 4) - Claculation of the outflow hydrograph at
:
[cfs]Now we have to use the curve 
to obtain 
from 
:
to obtain from :
[cfs] (interpolation)
- Calculation of the outflow hydrograph at
:
[cfs]
[cfs] (interpolation) - Computation continues for further time steps ⋯
The whole pocess can be formulated as follows:
- (1) Use the following discretized equation and compute the right hand side form data in the prevoisu time step:
- (2) Conduct a linear interpolation on storage-outflow curve for for
as follows:
- (3) Given
, compute the right hand side of Eq. 1 for the next time step as follows
- (4) Use the following discretized equation and compute the right hand side form data for the prevoisu time step and continue the iteration ...
Figure 21: The storage-outflow curve (left) and inflow and routed outflow hydrographs for the level-pool routing example.
4-2 Hydrologic River Routing
In order to calculate the routing for streams in a watershed, we can no longer assume that the water surface is level, since there is not negligible velocity in a stream and there are backwater effects. Therefore, we need to define a different storage-outflow equation, .
. The Muskingum method assumes that storage in a channel is the sum of a prism and a wedge storage.
- The prism storage is a linear reservoir defined as
where K is a constant representing approximately the travel time through the reach. - The wedge storage is defined as
, where X is a weighting factor ranging from 0 to 0.5 and thus,
This wedge storage represents a flood wave as it propagates through a stream channel as shown in Figure 22.
Figure 22: Schematics of wedge storage and flood wave propagation assumed in Muskingum method (Credit: The COMET Program)
We can discretize this storage equation as follows for time step jth
and thus:
Recall, we previously defined the change in storage as follows:
Combining equations 2 and 3, one can obtain the following flood routing equation:, known as the Muskingum stream routing formulation:
where 
.
. Typically, the streams is broken up into N subreaches when performing this calculation. The U.S. Army Corps of Engineers recommends the following condition to ensure solution stability when we discretize the river reach:
---------------------------------------------------------------
Example 7.4:The inflow hydrograph to a river reach is given below. Determine the outflow hydrograph of this reach, if 
[hr], 
and 
[hr]. The initial outflow is 85 
.
[hr], and [hr]. The initial outflow is 85 .
Figure 23: The results of the Muskingum routing.